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We present an approach for analyzing the dc current in voltage biased quantum 
superconducting junctions. By separating terms from different n-particle processes, 
we find that the n-particle current can be mapped on the problem of wave transport 
through a potential structure with n barriers. We discuss the relation between res- 
onances in such structures and the subgap structures in the current- voltage charac- 
teristics. At zero temperature we find, exactly, that only processes creating real exci- 
tations contribute to the current. Our results are valid for a general SXS-junction, 
where the X-region is an arbitrary non-superconducting region described by an 
energy-dependent transfer matrix. 
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1. Introduction 

Multiple Andreev reflection (MAR), first suggested by Klapwijk, Blonder, and Tinkham jj| for 
explaning the subharmonic gap structure (SGS) in SNS junctions, is presently accepted as a general 
mechanism of subgap current transport in superconducting junctions. During the last five years 
a considerable effort has gone into developing a consistent theory of MAR capable of including 
effects of quantum coherence and normal electron scattering by the junction. Different techniques 
have been used for calculating current-voltage characteristics, including various modifications of 
Keldysh formalism §, § |, § and the Landauer-Biittiker scattering method §, @, |, § |k| 0. 
The theory has been extended to resonant tunnel junctions (l^, [l4|, [l^] and junctions of d-wave 
superconductors Jl6t . The theory has successfully explained the SGS in transmissive planar junctions 
and atomic-size point contacts Jl8| . The good agreement between the theory and experimental 
data obtained for one-channel tunnel junctions without fitting parameters |ll| has provided a firm 
basis for further investigations - revealing the transport channels and transmissivity of individual 
channels in single-atom contacts [^0[ . 

Despite successful numerical calculations of subgap current- voltage characteristics for different 
types of junctions, the general understanding of MAR is still not sufficient, and interpretation of 
particular features of SGS, e.g. current peaks, presents difficulties. This especially concerns resonant 
MAR, where the energy dispersion of electron scattering phase shifts is important. 
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Fig. 1. Sketch of the junction geometry: The junction consists of an arbitrary normally conducting, single 
channel region (X), connected to two superconducting reservoirs (S). Two ballistic normal regions (N) of 
length e between the X and S regions are introduced for convenience, e is much smaller than the coherence 
length. 

In this paper, we present an approach where MAR is treated as a wave propagation problem 
in energy space (cf. ||). Using scattering theory formalism, we derive an adequate mapping of 
MAR onto a transmission problem through a ID wave guide with a built in multiple tunnel barrier 
structure. All parameters of this structure are uniquely determined by the junction characteristics 
(normal and Andreev scattering amplitudes) and by the applied voltage. In terms of such a mapping, 
the onsets and peaks of SGS are explained in terms of resonances in the transmission along the 
energy axis. The mapping allows us consistently to separate currents associated with n-particle 
transmission processes p2[ ] and to prove the cancellation of non-physical ground state currents, 
which is equivalent to the Pauli exclusion principle. 



We consider two superconducting reservoirs connected to a single normal conducting channel which 
may contain tunnel barriers (X-region) (see Fig. 1). The transport properties of this channel is 
described by its transfer matrix. Using an energy dependent transfer matrix we can model any 
effective single-particle potential structure, including the important scattering phase shifts in long 
and resonant junctions. In general the transfer matrix is also voltage dependent, and deriving this 
dependence for any given physical model is a problem in its own right. But having found this transfer 
matrix, we can determine the current when the contacting reservoirs are in the superconducting 
state. 

We place the origin of our coordinate system in the middle of the junction; the normal- 
superconducting (NS) interfaces are then located at —L/2 and L/2, where L is the length of the 
junction (Fig. 1). Due to the point-contact geometry of our junction, the superconducting pair 
potential can be considered steplike, A(x) = Ao0(|x| — L/2). 

We will construct scattering states in terms of the transfer matrix of the X-region. To this end we 
introduce two auxilliary normal regions, between the X-region and the superconducting electrodes. 
These regions are assumed to have the same material parameters as the electrodes, providing perfect 
NS-interfaces, and their length is much smaller than the coherence length. 

To calculate the scattering states we first make an ansatz in terms of a sum of plane- wave solutions 
to the Bogoliubov-de Gennes (BdG) equation in the left/right auxiliary normal regions (^ln/^rn) 
and in the left/right superconductors (^ls/^rs)- I n superconductors, the energy is commonly 
counted from the chemical potential which is motivated by the clcctron-holc symmetry This is 
not convenient in voltage biased junctions because the chemical potentials in the left ((J-ls) and 
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right (^rs) electrodes are different. A global reference of energy, which is particularly desirable 
in the case of an energy dependent transfer matrix, can be introduced by making different gauge 
transformations in the left and right superconducting electrodes. This procedure gives rise to the 
appearance of time-dependent phase factors in the wave functions of the superconducting electrodes 
implying inelastic scattering. To get symmetrical expressions for quasiparticles incoming from the 
left and from the right we choose (uls + Mhs)/2 as a global reference of energy, i.e. the reference 
energy is in the middle of the chemical potentials in the left and right superconductors. We thus get 
different time-dependences for electron and hole components in both superconductors, 
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In the ansatz in Eqs. E denotes the energy of the incoming particle, E n = E + neV, 

and the index a 6 {1,2,3,4} denotes the incoming quasiparticles {e^, h~*, e*~, h^} respectively, 
ufj = (u n ,v n ) and iv\ — (v n ,u n ) are vectors describing bulk electron- and hole- like quasiparticles 
respectively. 

By matching the wave functions across the NS intefaces we get the boundary conditions for 
the wave function in the normal region, including the source term. By defining coefficient vectors 
a n — (a n ,(3 n ) and j n = (j n ,S n ), and similar for the primed quantities in Eq. (^), we may write 
these boundary conditions in the following form, 



& n = U„7„_i + 5 n0 (ul ~ vl) f f_f 2 j° Q ) , l'n= Und^! + S n0 {ul - vl) ( ^Jjjl 1 . ( 5 ) 

where U n = diag(a n ,a~ r ) is the Andreev reflection matrix and a n is the amplitude of Andreev 
reflection, 

Vn \ (£„ - sign(£„ - A 2 ) /A \E\ > A 

a n = — = < v / , \ ' . (6) 

u » \ [E n - i VA 2 - El) /A \E\ < A 

In this derivation we have neglected the difference in wave vectors between the normal and super- 
conducting regions (quasiclassical approximation). The coefficients in Eqs. (|l|)-(^) are connected by 
the transfer matrix T{E) describing the X-region, 



Qn = T(£„ + eV/2)o! ri 



% = T(-(E n + eV/2))*/' n . 



(7) 
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Equations (|5|)~(]3)' together with the boundary condition of vanishing coefficients for \n\ — > oo, 
completely determine the scattering states. 

Checking equations (||)-(0) one may find that in any scattering state half of the coefficients in the 
ansatz will be zero. For example, for a S {1, 2} the nonzero coefficients will be a 2m & n d 72m-i ; while 
for a S {3, 4} a 2m _i and ^j 2m will be nonzero. To get rid of the redundant coefficients we define new 
coefficients c n ±. For a G {1, 2} the definition reads c 2m+ = a 2m , c 2m _ = 72m-i, c (2m+1)+ = % m+1 
and C(2 m +i)- = c\' 2rn . These coefficients are connected by transfer matrices T l n with even indices 

Tin = [T(£ 2 m + eF/2)]" 1 and odd indices T l 2m -x = T(-(-E 2 m-i + eV/2)). For the purpose of 
treating quasiparticles incoming from the left and right on an equal footing we introduce similar 
notations for a e {3,4}: c 2m+ = % m , c 2m _ = &' 2m _ x , c (2m+1 ) + = a 2m+1 and C( 2m+l) _ = j 2m and 

matrices T r 2m = T(-(E 2m + eV/2)) and TLj - [T(£ 2m _i + eV/2)] -1 . 

Using these definitions the equations (H)-® can i f° r quasiparticles incoming from both left and 
right, be written as: 

TT - I X I 2 2\ ( (fil(X + 5ia)/U0 \ - rrt - /Q\ 

c n + = U„c„„ + o n0 (u - v ) I _/g^ _|_ § 3 \/ Vq J ' c («+l)- = T„c n+ . (8) 

Note that even though these equations are formally identical for quasiparticles incoming from both 
left and right, the T„-matrices involved in the two cases are different. In Figure 2a we show the 
structure of Eqs. (^) for an electron-like quasiparticle incoming from the left. 

3. Mapping 

Now we are prepared to change focus from transport in real space to transport in energy space. 
The vectors c n ± have been defined so that the upper element is the coefficient for particles with 
positive k and the lower element is the coefficient for negative k. We notice that both electrons 
and holes with positive /c-value gain energy in passing the junction (for our choice of the sign 
of bias), while particles with negative k- value lose energy. Therefore the upper vector component 
describes upward motion in energy space, while the lower component describes downward motion. 
To emphasize this motion in energy space we introduce the following notation for the components 
of c n ±: 

where indicates upward motion in energy space and indicates downward motion (see Fig. 2a). 
With this notation, Eqs. (^) have obvious similarities to the problem of wave propagation through 
a one-dimensional multibarrier structure. Indeed the coefficient vectors are connected by simple 
matrix multiplication through source-free regions, n > m > or > n > m, 

U m+ iT m . (10) 

The M-matrices describe a ID multibarrier structure, where the T n -matrices provide tunnel bar- 
riers, and the U„-matrices introduce extra spacing, i.e. phase gained, between these barriers. 
Within the superconducting energy gap these M-matrices have the properties of ordinary real 
space Schrodinger transfer matrices, which conserve current: det(M) = 1, <7 z Mcr 2 = M'. In our 
case, the conserved quantity is, j n ± — \c] l± \ 2 — |c^ ± | 2 , which can be interpreted as the probability 
current flowing upwards in energy space. The probability current j n ± is related to electron and hole 
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Fig. 2. Schematic structure of a scattering state (a) and its probability currents (b). Figure (a) shows 
schematically the connection between different coefficients and transfer matrices for an electronlike quasi- 
particles incoming from the left (a — 1). Figure (b) shows schematically how the probability current from 
the incoming quasiparticle enters the normal region and divides into upgoing (jo+) arid downgoing (jo-) 
currents. This probability current then leaks out due to incomplete Andreev reflections outside the gap 
and/f). 



currents in the original real space problem. For example, for quasiparticles incoming from the left 
i(2n+i)+ = i(2n+2)- is probability current of electrons, while —jm-v = — j(2n+i)- is probability cur- 
rent of holes (the equalities are provided by the current conserving properties of the normal transfer 
matrices T). The quasiparticle probability current 1%, which is conserved by the BdG equation, 
consists of the sum of electron and hole probability currents, 1% — j(2«+i)+ + ( — j(2n+i)-)- The 
probability current I? is obviously equal to zero inside the energy gap because of the conservation 
of the current j n . 

The M- matrices are associated with scattering matrices S, 

4- ) nm Ui- ) ' nm { r- m ) W 

where the reflection amplitudes r + = —M2i/M22,r~ = M12/M22 and transmission amplitude t = 
I/M22 satisfy the current conservation condition |r| 2 + \t\ 2 = 1 within the superconducting energy 
gap. Outside the gap this condition changes into |r^„J 2 + Knm| 2 Ilp=m+i l a pl~ 2 — 1j where the 
equality holds only for fully transparent junctions. Since the probability of Andreev reflection is 
smaller than unity outside the gap, |a ra | 2 < 1 for \E n \ > A, the current j n is not conserved, which 
is related to leakage outside the normal region due to incomplete Andreev reflection. This leakage 
is characterized by the the probability current I? 7^ 0. The sign of I? in the above definition was 
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Fig. 3. Schematic structure of the mapping of the scattering state equations on a leaky one-dimensional 
Schrodinger potential. The T„-matrices introduce tunnel barriers (black boxes). Inside the energy gap, 
\E n \ < A, the Un-matrices correspond to extra spacing (white box), while outside the gap they introduce 
leaking of particles out of the channel. The sources (5 a i + Soi and 8 a i + 5 a i) introduce particles into the 
potential structure and the amplitude for particles going away from the source are G„o and GW_„\, n > 0. 
The particles are reflected from ±oo with reflection amplitudes r„±. 



chosen so that probability current flowing away from the junction is positive, and since the only 
incoming current is from the source, I? > for n ^ 0. We will see in the next section that the dc 
current can be expressed in terms of this "leakage" . 

Particles are introduced into this "leaking" potential structure by the source. It follows from 
Eq. (||) that the injection of quasiparticles with positive k (a £ {1,4}) introduces upgoing particles 
directly above Uo, while quasiparticles with negative k {a € {2,3}) inject downgoing particles just 
below Uq (see Fig. ||). 



4. Formal solution 

To solve Eqs. (||) one has to find two independent homogenous solutions that go to zero for n — ±oo 
respectively, and then to match them at n = 0. The boundary conditions, 

lim M nm c m+ = (m > 0), lim (M„ m )- 1 c„_ =0 (n < 0), (12) 



m — > — oo 



fix the ratios cl + /cl + = r n+ = lim m _ >00 r+ n for n > and c„_/c„_ = r„_ = lim ln _ ) ._ 00 r nm for 
n < 0. The quantities r n ± have the meaning of reflection amplitudes from ±oo for upgoing and 
downgoing particles respectively. 

We can now express the scattering state coefficients in terms of the matrices S nm and the reflection 
coefficients r n ±. We choose to give expressions for c„_ for n > and c n+ for n < 0, knowing that 
Eq. (^) connects c„+ to c„_. 

• [S la + S 2a a Q r l Q _] ■ G' rM ■ ( 2 \ ) (n > 0), (13) 



Uo 

_ "li " ' li I:" . / .-./ / 

Cn+ — 

U 



^ • [6 2 . + S lt7 a r 0+ ] ■ G l Qn ■ ( ) (n < 0). (14) 



The expressions for a £ {3,4} is found from Eq. (13) by the substitutions I — * r, 5\ a — ► 8± a and 



82a ^3<r- The expression for the quantity G in the above equation reads 



Superlattices and Micro structures, Vol. ??, No. ??, 1999 



7 



Gnm = n _ 2 + vi _ 2 ~ - ^_ & — ~ ( n > TO )- ( 15 ) 

G nm is the dressed propagator, through source free regions of the multibarrier structure from point m 
to point n (see Fig. 0). The difference between the dressed propagator G nm and the bare transmission 



amplitude t nm is due to reflections from the region outside [E m ,E n ]. Note that G nm in Eq. (15) is 
expressed through r n+ and r m _ (using the fact that the reflection amplitudes from the same sign 
of infinity are related by Eq. (0)). This form is vital for the cancellation theorem below (Eq. (p2|)), 

5. dc Current 

The total dc current through the junction has the form, 

I**=\\ dE [ E 1 = f{E)Y,f(E), (16) 
hJ\ E \>A VE 2 - A 2 ^ 

where the function f{E) denotes the equilibrium population of the scattering states, originating 
from bulk electrodes. The current density j a (E) has the form, 

oo 

f{E)= J2 (\<(E)\ 2 -\(3Z(E)\ 2 + h:(E)\ 2 -K(E)\ 2 ). (17) 

n— — oo 

We can rewrite this in terms of the current flowing upwards in energy-space, j n ±, 



oo 

J%E) = X>n+(£) + E 3n-(E). (18) 

n— n— — oo 

We note that for positive n the current j n ± is positive, while for negative n it is negative, and 
therefore the two terms in the above equations have opposite signs. Using the boundary condition 
of vanishing coefficients for large |n|, we can rewrite the summation over j n ± in Eq. ( |l8| ) as a 
summation over probability currents I^(E), 

f(E) = J2nI^{E). (19) 

njtO 

We have now managed to separate the dc current into a weighted sum over probability currents 
transmitted from Eq to E n , where the weight n, is the number of times the probability current 
passes the junction. Since each Andreev reflection at energies within the superconducting gap is 
associated with the transfer of charge 2e, Eq. ( |l9| ) presents the total charge current as a sum over 
n-particle transmission processes p2[ . 

The total probability current that flows between E and E n , including the superconducting density 
of states, is 



Z " (g)= ^- A ^ ™ = £ (l"l«o| 2 )(l-|a„| 2 )|G« | 2 (l + K< + | 2 )(l + |a r Q _| 2 ). (20) 

f a£{l,r} 

This current has no divergences, moreover by using the conservation law for probability current one 
finds that I%(E) < 4. The form of the current in Eq. ( pp| ) can be interpreted in the following way. 
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The current is proportional to the probability of transmission from E to E n , and the terms (1 — |ao| 2 ) 
and (1 — |an| 2 ) are the probabilities to enter and leave the normal region respectively. Transmission 
from E to E n consists of four alternative paths: direct transmission with the probability |G n o| 2 j 
transmission with excursion to and reflection from either minus infinity or plus infinity, which yields 
an additional factor |aoro_| 2 or |a„r n+ | 2 respectively, and finally transmission with excursions to 
and reflections from both minus infinity and plus infinity. 

The probability current I%(E) in Eq. (|2^) obeys the important equation 



I p _ n (E)=I%(E-neV), (21) 

which implies a detailed balance between probability currents flowing upwards and downwards in 
energy space respectively. This balance is straightforward to prove using the equality between S nm 
for quasiparticles incoming at E + qeV and St n +q\r m + q \ for quasiparticles incoming at E. Using this 
detailed-balance equation one is able to prove the exact cancellation, at T — 0, of currents which 
do not pass the energy gap, 



\-neV r-A 

P n (E)dE- I p _ n (E)dE = (n>0). (22) 

This cancellation theorem proves the consistency of the single particle scattering approach to su- 
perconducting junctions. It is well known that within the Landauer-Biittiker scattering approach 
to elastic tunneling, e.g. in normal junctions |^3|, the Pauli exclusion principle is automatically 
valid due to the exact cancellation of currents that do not produce real excitations. The cancellation 
theorem Eq. ( |22|) extends this property to inelastic scattering (MAR) in superconducting junctions. 

Using Eq. (|22|), and also the electron-hole symmetry of the BdG equation (J^ j CT (~ E) = 
— j a (E)), we arrive at the final formula for the dc current at finite temperature, 

(/.-A p-A-neV \ 

Q(neV - 2A) / dE[2f(E) - 1]I*(E) + 2 / dE[f(E) - f(E n )]I*(E) . 

JA-neV J-co j 

(23) 

The obtained formula for the dc current in Eq. ( p3[ ) has only positive terms and consists of two parts. 
The first part is related to the creation of real excitations, due to quasiparticles traversing the energy 
gap, and it is responsible for the SGS. This current exists at zero temperature and decreases with 
increasing temperature. The second part is a smooth background current from thermal excitations, 
and is exponentially small at low temperatures. 



6. Subgap structure 

Let us analyze SGS using Eqs. (jl^), ( pcf ) and (p3|). The magnitude of the n-particle current I%(E) 
is proportional to the transmittivity \t n o\ 2 of the n-barrier tunnel structure, which is estimated as 
|^no| 2 ^ D n in the absence of resonances, where D is the transmissivity of the normal junction. 
Therefore, the current in Eq. ( p3| ) consists of a step-like structure with the onsets at eV = 2A/n. 
This step- like structure is particularly pronounced in junctions with low transmittivity D <C 1 0, 
and it is washed out in transparent junctions with D w 1. However, even in the latter case, the 
current decreases exponentially at low voltages, Idc oc D 2A / eV , as soon as D ^ 1 B. 

This simple step-like structure is complicated by resonances in I^{E). One may distinguish three 
different sources of resonant behaviour. First there may be normal electron resonances in the transfer 
matrices T n 
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Secondly, even in the absence of normal electron resonances there are specific superconducting 
resonances in the bare transmission amplitudes i„o due to electron-hole dephasing during transmis- 
sion through long junctions and during Andreev reflections given by matrices U. In the mapping, 
dephasing is a source of phase gained between the barriers in Fig. || In short constrictions, L = 0, 
the resonant condition reads arg{a^) = rrnr, where m is integer, and it is fulfilled near the gap edges, 
Ek ~ ±A. (For energies outside the gap the phase gained is constant, ±7T, but for |-E/A| — 1 > D 2 
the resonance vanishes because of the leakage.) The resonance effectively removes two barriers so 
the bare transmission probability is increased by the factor D~ 2 . It is easy to see that at voltages 
eV = 2A/n it is possible to have simultaneously two resonances with corresponding indices related 
as k' — k + n. In this case the transmission probability \t n o\ 2 is enhanced by the factor D~ 3 if the 
resonances are next to each other (n — 1), otherwise the enhancement factor is D~ 4 . 

Finally there are the resonances due to reflection from ±oo in the denominator of Eq. (fl5|), 
boundary resonances. The conditions for resonance are (1 — ao r o- r no) = an< ^ (1 ~~ a n r n+ r no) = 0> 
which is possible for E = — A or E n = A respectively. These resonances are only important when 
they overlap with resonances in the bare transmission amplitude. An interplay of the resonances in 
the bare transmission amplitude t n0 and the boundary resonances produces current peaks in the 
subharmonic gap structure (Fig. ^ . 

As an example of this interplay let us discuss the 4-particle current in a short junction which 
possesses all typical features. The 4-particle current (I4) (see Fig. ||) has an onset at 2A/eV = 4. 
In the voltage region 3 < 2A/eV < 4 there are no resonances in £40 and therefore the magnitude of 
the current is I4 cx D 4 . For voltages 2A/eV < 3 a single resonance in t4o becomes possible, which 
gives an enhancement of the current to I4 cx D 4 \ln(D)\. Close to 2A/eV = 3 this resonance overlaps 
with a boundary resonance which results in a current peak with the magnitude {Ii) ma x ^ D 3 . At 
2A/eV = 2 there is a double resonance, E\ — —A and E 3 — A, which results in a current peak 
of order (l4) max cx D 2 . The voltage 2A/eV — 1 is a special case because the two resonances are 
next to each other so the enhancement of the transmission probability is only D~ 3 , resulting in a 
current peak (l4,) max cx D 3 . At this voltage the 4-particle current is weakened by the leaky Andreev 
reflection outside the gap. 



7. Conclusions 

In conclusion, we have presented an approach, for analyzing the SGS in quantum superconducting 
junctions, where MAR is treated as a transport problem in energy space. Such an approach is 
natural for the inelastic scattering problem represented by MAR. Using scattering theory formalism, 
we have derived a mapping of MAR onto a problem of transmission through a ID wave guide with 
a built in multiple tunnel barrier structure. There is no current conservation in this transmission 
problem due to a leakage outside the wave guide. The charge current in the original superconducting 
junction is expressed through this leakage, which is identical to the probability current of outgoing 
quasiparticles in the superconductors. In terms of such a mapping, the SGS is explained in terms of 
resonances in the transmission along the energy axis. Three different types of resonances have been 
found which are important for the interpretation of SGS: (i) superconducting resonances induced by 
electron-hole dephasing and the Andreev reflections, (ii) boundary resonances induced by reflections 
from the boundaries of the wave guide, and (iii) normal electron (Breit-Wigner) resonances, which 
may exist in addition to former two types of the resonances, which are always present in MAR. 
The mapping allowed us to separate currents associated with n-particle transmission processes and 

f The more complex cases of junctions with Breit-Wigner resonances and d-wave junctions are discussed in 
Refs. |15| and [E5[ respectively. 
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Fig. 4. The SGS of a SIS junction with transparency D = 0.01 divided into different n-particle currents, for 
n 6 {1,6}. The inset shows the SGS of a SIS junction with transparency D = 0.1 close to eV = A and the 
contributions from 2-, 3- and 4-particle currents. 



to prove the cancellation of non-physical ground state currents, which is equivalent to the Pauli 
exclusion principle. 
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